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Electrophysiological studies have established that the permeation of Ba^^ ions through the KcsA K^-channel is im- 
peded by the presence of ions in the external solution, while no effect is observed for external Na^ ions. This 
Ba^* "lock-in" effect suggests that at least one of the external binding sites of the KcsA channel is thermodynami- 
cally selective for K^. We used molecular dynamics simulations to interpret these lock-in experiments in the context 
of the crystallographic structure of KcsA. Assuming that the Ba^^ is bound in site Sg in the dominant blocked state, 
we examine the conditions that could impede its translocation and cause the observed "lock-in" effect. Although 
the binding of a K^ ion to site Si when site S2 is occupied by Ba^^ is prohibitively high in energy (>10 kcal/mol), 
binding to site So appears to be more plausible (AG > 4 kcal/ mol) . The 2D potential of mean force (PMF) for the 
simultaneous translocation of Ba^^ from site S2 to site Si and of a K^ ion on the extracellular side shows a barrier 
that is consistent with the concept of external lock-in. The barrier opposing the movement of Ba^^ is very high 
when a cation is in site So, and considerably smaller when the site is unoccupied. Furthermore, free energy pertur- 
bation calculations show that site So is selective for K^ by L8 kcal/mol when S2 is occupied by Ba^*. fiowever, the 
same site So is nonselective when site S2 is occupied by K^, which shows that the presence of Ba affects the selectiv- 
ity of the pore. A theoretical framework within classical rate theory is presented to incorporate the concentration 
dependence of the external ions on the lock-in effect. 



INTRODUCTION 

channels are a broad family of membrane proteins 
that are present in almost every cell. Their high selectiv- 
ity is one of the most remarkable aspects of cellular 
physiology. While remaining highly selective for ions 
over Na^ ions, these channels allow ions to cross a 
cellular membrane through a passive mechanism at a 
rate that nearly matches bulk diffusion. Conduction oc- 
curs through a transmembrane domain with tetrameric 
architecture. The subunits of this tetrameter contain a 
highly conserved sequence of amino acids, TTVGYGD, 
which is essential for selectivity (Heginbotham et al., 
1994). The high-resolution crystal structure of KcsA, a 
pH-gated K^ channel from the bacteria Streptomyces livi- 
dans, revealed that the selectivity filter provides a nar- 
row pore for single-file ion conduction comprised of a 
series of ion binding sites where the ion is coordinated 
by backbone carbonyl oxygens (Doyle et al., 1998) . This 
structure is illustrated in Fig. 1. The robust selectivity 
achieved by such a structurally simple motif is even 
more intriguing considering the sizable structural fluc- 
tuations that proteins undergo at physiological temper- 
ature (Allen et al., 2004). This is especially remarkable 
given that the sizes of K* and Na* ions are very similar; 
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the ionic radius of K^ is only ^^0.4 A larger than Na* 
(Shannon, 1976) . 

Under physiological conditions, the open-conductive 
state of the channel must allow K^ ions to move out- 
wards from the intracellular solution to the extracellu- 
lar solution without allowing Na^ ions to enter from the 
extracellular solution. In effect, this means that the nar- 
row filter must select for K* ions in its open-conductive 
conformation, without undergoing a long-timescale con- 
formational change. 

Since the report of the high-resolution x-ray struc- 
ture, KcsA has become a prototypical model for experi- 
mental as well as computational studies of the structure 
and activity of K* channels. A series of perspectives on 
ion channel selectivity was recently published in this 
journal (Andersen, 2011). Many computational studies 
have attempted to explain the structure, dynamics, and 
selectivity of the filter (Sansom et al., 2002; Shrivastava 
et al., 2002; Bemeche and Roux, 2003; Domene et al., 
2008; Domene and Furini, 2009; Furini and Domene, 
2009) . In particular, there have been many efforts to 
elucidate the origin of selective ion conduction using 
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molecular dynamics simulations (Allen et al., 2000; 
Aqvist and Luzhkov, 2000; Shrivastava et al., 2002; Noskov 
et al., 2004; Roux, 2005; Egwolf and Roux, 2010; Jensen 
et al., 2010; Kim and Allen, 2011). These studies have 
led to a range of interpretations of ion channel selectiv- 
ity, with Roux and coworkers concluding that selec- 
tivity is enforced by a higher free energy for Na^ ions 
binding in site S2, whereas Kim and Allen (2011) con- 
cluded that no binding site is significandy selective for 
K*. In part, the myriad of theories on the origin of the 
K^ channel selectivity result from difficulty in relating 
these simulations to the experimentally observed activ- 
ity of these channels. 

One of the most direct measures of K^ selectivity is 
provided by blockade experiments in the presence of 
Ba'* ions (Hagiwara et al., 1978; Armstrong and Taylor, 
1980; Neyton and Miller, 1988; Vergara et al., 1999; 
Piasta et al., 2011). As the radius of Ba^* is very similar to 
that of K"^ (1.38 A and 1.35 A, respectively) , both of these 
cations can fit inside the binding sites and permeate 
through the selectivity filter of a K^ channel. However, 
the divalent Ba^* ions bind more strongly to the selectiv- 
ity filter than K^ ions and, thus, permeate through the 
channel at a much slower rate. The binding of Ba^^ to 
the filter prevents the passage of other ions, resulting in 
blockades of the channel current on the millisecond 
timescale that are observed in single-channel record- 
ings. Recendy, such Ba^^ blockade experiments were 
performed for the first time with the KcsA channel (Piasta 
etal., 2011). Such experiments are expected to be particu- 
larly informative in this case because of the availability 
of high-resolution x-ray structure of this channel. When 
a 30— 60-pM concentration of Ba^^ ions was introduced 
in the "internal" solution, the single channel current 

extracellular 




cavity 

Figure 1. The structure of the selectivity filter of the KcsA 
K* channel (PDB accession no. 1K4C). The configuration where 
K* ions (green spheres) occupy sites So, Sq, and S4 is depicted. 
Only backbone atoms and coordinating threonine hydroxyls are 
depicted. Only one opposing pair of the tetramer subunits is 
shown for clarity. The cavity is formed by the gating helices in the 
closed state and is filled with water molecules. The extracellular 
solution is immediately above the selectivity filter. 



underwent frequent interruptions due to Ba^^ binding 
to the selectivity filter. The block is relieved when the 
Ba^^ ion exits the selectivity filter, either by returning to 
the intracellular side or proceeding to the external solu- 
tion. Block times were typically of the order of 1 ms in the 
absence of external K^, but increased to '^100—1,000 ms 
with an external K^ concentration of 5 mM. Even '^5-]aM 
external concentrations of K^ increase the duration of 
the Ba'* blockade lifetimes, which suggests the existence 
of an outer binding site with a high affinity for K^. In 
contrast, experiments with Na* ions in the external solu- 
tion showed no lock-in effect at all at concentrations up 
to 0.2 M. This dependence of the block time on the con- 
centration of extracellular ion is interpreted as evidence 
of the so-called K^ lock-in effect, where the binding of 
K^ in a site external to the Ba^^ prolongs the blocking 
time by impeding its translocation forward toward the 
extracellular side. Because the blocking times are very 
long, these experiments can be interpreted by assuming 
that binding of a cation to a site located external relative 
to the bound Ba^^ is in near thermodynamic equilibrium 
with the extracellular bulk solution. Using a kinetic model 
based on this view, it was deduced that the external lock- 
in site must be thermodynamically selective for K^ over 
Na^ by at least 7 kcal/mol. 

In many respects, these Ba^^ blockade experiments 
are highly complementary to detailed free energy com- 
putations using molecular dynamics simulations. As the 
permeation of Ba^* occurs on a time scale that is many 
orders of magnitude larger than the rate of Na^ or K* 
association and dissociation from the filter, the alkali 
ions binding in the external sites of the filter can be 
safely understood as being in equilibrium with the ex- 
ternal solution. In this way, the selectivity of the channel 
does not need to be discussed in terms of the nonequi- 
librium diffusive process of conduction, but rather it 
is framed in terms of the equilibrium thermodynamic 
binding affinity of the external sites to an alkali ion and 
the quasi-equilibrium process of rare, activated transi- 
tions of Ba^* between sites. The Ba^* blockade experi- 
ments form the basis of a thermodynamic view of ion 
selectivity, where differences in the binding Etffinity of K* 
or Na^ for specific sites along the narrow pore is the 
mechanism by which these channels achieve selective 
ion permeation. In this paper, we report umbrella sam- 
pling and free energy perturbation (FEP) molecular 
dynamics simulations designed to model the Ba^^ block- 
ade experiments of Piasta et al. (2011) in an attempt to 
reconcile the experimental observations with the re- 
sults of simulation. 



MATERIALS AND METHODS 

Models and computational parameters 

All simulations reported here were performed using the pro- 
gram CHARMM, version c36b2 (Brooks et al., 2009). The atomic 
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simulation systems were based on tlie systems that were used in 
previovis simulations by our group (Berneche and Roux, 2000; 
Egwolf and Roux, 2010). The tetrameric crystallographic struc- 
ture of the KcsA channel (Protein Data Bank [PDB] accession no. 
1K4C) reported by Zhou et al. (2001) is embedded in a 70 A x 
70 A bilayer comprised of 112 dipalmitoyl-phosphatidylcholine 
(POPC) lipids. The pore axis of the channel was aligned along 
the z axis of the simulation box, normal to the bilayer This pro- 
tein-membrane layer was solvated by 6,778 water molecules, 
forming a bulk liquid solution around the membrane. CI" and 
K* ions were introduced into the bulk solution to neutralize the 
charge of the system and establish an ionic concentration equiva- 
lent to a 150-mM solution of KCl. As the number and valency of 
the ions in the filter were adjusted, the net charge of the system 
was maintained at neutrality by adjusting the number of CI" and 
K* ions in solution. 

The electrophysiological experiments of Piasta et al. (2011) 
used the E71A KcsA mutant to prevent their measurements of the 
Ba^* blockade events from being obscured by C-type inactivation, 
which does not occur in the E71A mutant. This type of inactiva- 
tion occurs on a much longer timescale than the length of our 
simulations, and no conformational change corresponding to 
this type of inactivation was observed in our simulations. As the 
E71A mutant otherwise shows the same conductivity as the wild- 
type channel modeled in this report, our simulations of Ba^* 
blockades of the wild-type structure should generally be applica- 
ble to this mutant as well. 

All simulations performed here made use of the CHARMM 
PARAM22 protein force field (MacKerell et al., 1998) with the 
backbone dihedral Cmap potential correction (MacKerell et al., 
2004). The CHARMM PARAM27 force field force was used for 
the lipid parameters (Feller et al., 1997). Optimized Lennard- 
Jones parameters for Na"^ and K* were used (Noskov et al., 2004). 
The parameters for Ba^+ (£„„ = -0.150 kcal/mol, R„„ = 1.849 A) 
were determined by adjusting the Lennard-Jones radius of Ba^* 
to match the experimental hydration free energy. These param- 
eters were tested by comparing the computed radial distribution 
function (RDF) calculated using these parameters to a QM/MM 
molecular dynamics simulation (see Rowley and Roux [2012] 
and Riahi et al. [2013] for the details of this type of simulation). 
The RDF maxima of these two methods were identical within 
0.1 A, indicating that the Ba^* parameters are reasonable. All 
other Lennard-Jones interactions were calculating using the com- 
bination rule, with the exception of the Ba^* — O(carbonyl) in- 
teractions (E„in [Ba^* - O] = -0.134 kcal/mol, R„i„ [Ba^+ - O] = 
3.36 A), which were adjusted to reproduce the relative R1MP2/ 
def2-TZVP hexacoordinate ion-ligand binding energies of N-methyl- 
acetamide and water Water molecules were described using the 
T1P3P model (Jorgensen et al., 1983). Bonds containing hydrogen 
atoms were constrained using the SHAKE algorithm (Ryckaert 
et al., 1977). The electrostatic interactions were computed with 
the particle mesh Ewald (PME) method, with a72Ax72Ax8lA 
grid (roughly 1 grid point per angstrom; Essmann et al., 1995). 
The systems were simulated with a time step of 2 fs. The tem- 
perature and pressure of the system was regulated by the CPTA 
method (Feller et al., 1995, 1997), where the surface area of the 
membrane in the xy plane was kept constant while the length 
of the unit cell along the z axis was allowed to vary to preserve a 
constant pressure. 

Umbrella sampling simulations 

The umbrella sampling simulations of the potential of mean 
force (PMF) were performed by applying planar harmonic bias- 
ing positions to the Z coordinates of the translocating ions. The 
details of this method are presented in a paper by Berneche and 
Roux (2001). The 2D PMFs were calculated by an umbrella sam- 
pling simulation where both the Z coordinate of the Ba^* ion ( Zsa) 



and the Z coordinate of the lock-in ion (Zl) were simultaneous 
restrained. An initial grid of ion positions along the Zbj axis was 
defined based on the locations of the ion binding sites in the 
1K4C crystal structure. For translocation of Ba^* from site Sg to Si, 
this grid covered the range Zga = 0-5.5 A in increments of 0.5 A. 
The range of 2^ positions was initially assumed to extend from 
Zbb = 3.5 A to 14 A. To generate stable structures from these artifi- 
cially constructed .systems, energy minimizations were performed 
where the positions of the ions, bulk water molecules, lipids, and 
most of the proteins were fixed. The oi-carbon of the selectivity 
filter residues was restrained by a 5-kcal A"" harmonic force and 
the water molecules within the channel were unrestrained. These 
systems were minimized by the adaptive Newton-Raphson algo- 
rithm for 200 steps. The constraints were then removed and 
300-ps-long molecular dynamics simulations were performed on 
these .systems to equilibrate the structure with the Ba^* and lock-in 
ions restrained to their assigned positions. A spring constant of 
10 kcal moP' A"^ was used for the planar harmonic restraint on 
the lock-in ion, while a stronger spring constant of 50 kcal moP' A"^ 
was used to restrain the position of the Ba^* ion. The lock-in ion 
was restrained to remain inside a cylinder with a 3-A radius 
around an axis running along the center of mass of the filter 
using a flat-bottom half harmonic restraint with a force constant 
of 10 kcal A"l 

The ID PMFs of K* binding to the external sites of the filter 
when site Sj was occupied by K* or Ba^* presented in Fig. 2 were 
calculated by integrating over the regions of the corresponding 
2D PMFs that correspond to Ba^* being bound in site S? (0 A < 

^a<2A). 

2A 
OA 



Hamiltonian-replica exchange MD 

To improve the rate of convergence and ergodicity of our con- 
figurational sampling of the umbrella sampling windows used 
to calculate the PMFs, we used Hamiltonian replica exchange 
molecvilar dynamics (H-REMD) to allow exchanges between win- 
dows. In this method, attempts are made to exchange the config- 
urations of neighboring systems on the PMF, with the acceptance 
of these attempts governed by the Metropolis algorithm (Jiang 
and Roux, 2010). This type of strategy has been demonstrated to 
significantly improve the convergence of distributions sampled 
using molecular dynamics simulations. 

The 2D PMFs calculated in this study required an REMD 
scheme that is more elaborate than usual because umbrella-sam- 
pling windows are not distributed in a Cartesian grid, where the 
nearest neighbors provide a natural list of partners for exchanges. 
Additionally, because of the steepness of the free energy surface 
along the Zb^ coordinate, the restraining forces on the Ba^* ion 
are large (50 kcal moF' A"'), leading to low exchange rates be- 
tween these windows. Although reducing the spacing of the 
windows could improve the exchange rate, this would greatly in- 
crease the number of windows. To resolve these problems, we used 
a novel strategy we term "woven" H-REMD. In this method, a ID 
Hamiltonian exchange sequence is constructed from the 2D grid 
of windows by all windows with given values of in a ID sequence. 
At the maximum or minimum value of Zk for a given value of Zg,, 
exchanges are also attempted with a neighboring replica with an 
adjacent value of Zbb. This strategy allows irregularly shaped um- 
brella sampling grids to be sampled within the existing Hamiltonian 
replica exchange molecular dynamics method implemented in 
CHARMM. To limit the number of replicas, the windows were 
separated by increments of 0.5 A, which still allowed for signifi- 
cant overlap and an average exchange acceptance probability of 
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20%. In total, the 2D PMFs for the permeation of Ba^* from site S2 
to site Si reqviired 162 replicas. 

The unbiased 2D PMFs were computed from the umbrella 
sampling simulations by the weighed histogram analysis method 
(WHAM; Roux, 1995). The histogram bin width was 0.1 A. The 
WHAM equations were iterated until convergence was achieved 
within a tolerance of 10"^ for all windows. 

Definition of 1 D translocation PMF 

Let t/(X) represent the potential energy as a function of all atomic 
coordinates X in the system, and n represent the number of ions 
bound to the outer sites of the pore. Let //„(X) be an indicator func- 
tion equal to 1 when the system is in state n, and zero otherwise. If 
we have a complete set of indicator functions, then ~ ^ by nor- 
malization. For the binding of a monovalent cation to the outer sites 
of the selectivity filter, we have only two states, with or without a 
bound cation. For the average of a property A(X), 



(A) = ZR) 



(2) 



where P„ = (//„) is the probability the state n, and (A)(„) repre- 
sents the conditional average, 
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Let / represent the position of the barium ion along the chan- 
nel axis. For a Z-constrained average of a property A(X), the ex- 
pression is a little more complicated. 
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Here, X' is the integration variable and Z is the fixed 
constraint. 
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= Z(4,z,„)^»(2). 

n 

The mean force along the Z coordinate is then, 



For the Ba^* lock-in effect, we consider two states: n = 0, where 
there is no lock-in ion, and n = 1, where there is, 

(^')(z)=(Az,»^o)^o(2) + (F)(,„^„Pi(Z) 



where K,giZ) is the equilibrium binding constant of the outer cat- 
ion when Ba^* is held fixed at Z. A delta function is used to restrict 
the configurational integrals to the holo and apo states. The equi- 
librium binding constant of the outer ion should be written as. 
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This expression is valid for any arbitrary position Z. Hence, 
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\dZ' 
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FEP 

The free energy of substitution of K* and Na* (AAGK-^Na) was cal- 
culated using FEP molecvilar dynamics (FEP/MD) as the differ- 
ence between the free energy to replace with Na* in a binding 
site and the free energy for the same transformation in bulk water 
(AAGK-^Na = AGK_^Na.site - A GK->Na,buik) • The new potential energy 
function is defined as a linear combination of the K* ( L^k) and Na* 
containing ( t5.ja) potentials according to the formula 

U{X) = {l-X)U^ + XU^^. 

The free energy difference of these two systems was calculated 
using the relation 



AG, 



The FEP was performed using a total of 1 1 values of X in incre- 
ments of 0.1 between 0 and 1. Each window was equilibrated 
for 200 ps before a production sampling period of 500 ps. The 
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integration of these windows to determine the free energy differ- 
ence and the sampling error was performed using the multistate 
Bennett acceptance ratio (MBAR) method (Shirts and Chodera, 
2008) . The lock-in ion was restrained to remain inside site So by 
flat-bottomed planar harmonic restraints, which have zero force 
when the / coordinate of the ion is inside the site (defined by the 
carbonyl ligands at the top and bottom of the site) but experience 



a strong quadratic restraining force (A„ 
the ion moves outside the site. 



100 kcal mor' A"-^) if 



Calculation of the transmission coefficient 

We estimated the rate of Ba^* permeation from site S2 to site Si 
when there is no lock-in ion using Grote-Hynes rate theory 
(Grote and Hynes, 1980). Eyring transition state theory neglects 
dissipative effects that diminish the reaction rate. This can be cor- 
rected for by calculating the reaction transmission coefficient, k, 

The transmission coefficient can be computed in a straight- 
forward fashion using Grote-Hynes theory. 



where o), is the frequency corresponding to the cuwature of the 
PMF at the top of the barrier, i.e., = |[w"(Z*||^m, and s is the 
"reactive" frequency corresponding to oscillations across the barrier. 
This reactive frequency can be determined from the equation 

.2 2 >^M(s) 



where M{^s) is the Laplace transform of time-dependent friction 
at the top of the barrier. 



M{s) = \dte^"M{t). 



M{^s) can be determined by a technique proposed by Straub 
et al. (1988) based on analysis in terms of the generalized Langevin 
equation for a harmonic oscillator 



c{t)- 



8Z^ 



-\dt' C{t')- \dt' M{t-t')C{t'). (3) 



Where Z is the velocity of the ion along the reaction coordinate 
and 6Zis the displacement of the system from the transition state 
(SZ= Z — Zfs), and C(t) is the velocity autocorrelation function. 
The Laplace transform of Eq. 3 gives 



KT C{s) 



M{s)C{s). (4) 



By substituting C(0) = ^Z^^ and m = kgT / \^'].^^, Eq. 4 can be 
simplified to, 



M(s) 



5zn 



(5) 



Using Eq. .5, M{s) / m can be calculated using a molecular 
dynamics simulation that is tightly restrained to the top of the bar- 
rier. (-Z^^^ and l^Z^^^ are calculated from ensemble averages 



of this trEijectory. £(5) is calculated by a numerical Laplace trans- 
form of the velocity autocorrelation function from the molecular 
dynamics simulation. 

The well frequency was determined by approximating the no 
lock-in PMF in Fig. 6 as a parabola with a curvature of W" (Z*) = 
25 kcal moL' A"^. To determine {z^}ts, (S^}ts, and the veloc- 
ity autocorrelation function (C(i)), a .500-ps molecular dynamics 
simulation was performed where the Ba^* ion was restrained at 
the top of the barrier with a 50 kcal mol"' A"^ harmonic restraint. 
These data were used to calculate M(i) numerically over a range 
of values of i. The reactive frequency, s, was determined through 
a numerical solution of Eq. 4 using the value of 0)+ from the PMF 
and these values of M(i), yielding a reactive frequency of s = 3 ps~^ 
and a transmission coefficient of k = 0.25. 



RESULTS AND DISCUSSION 

PMF of lock-in ion binding 

Piasta et al. (201 1) used the available structural and elec- 
trophysiological data for KcsA to assign the Ba"^ binding 
and lock-in sites. Analysis of the K^-free block time dis- 
tribution was consistent with two distinct Ba"* binding 
sites within the filter. Because the electronic density as- 
sociated with Ba^^ was only observed in sites S2 and S4 in 
the crystal structure of KcsA soaked with Ba^^ (Lockless 
et al., 2007), it was assumed that sites 89 and S4 are the 
sites that are predominantly occupied by Ba^^ during 
the blocking events. Based on the analysis of electro- 
physiological data, it was inferred that the long blocks 
occur predominantly when the Ba^* is bound to site S2, 
implying that the external cation must be binding to 
either site So or site Si. 

Piasta et al. (2011) argued that the high selectivity dis- 
played in the experiments is inconsistent with So being 
the lock-in site, and postulated that site Si must be the 
actual lock-in site. One may note that this leads to a 
multi-ion binding arrangement that differs from the 
generally accepted view, in that two ions are simultane- 
ously bound to adjacent sites (S2 and Si). X-ray crystal- 
lography of ion-bound channels (Zhou and MacKinnon, 
2003) , streaming potential measurements (Alcayaga et al., 
1989) , and molecular dynamics simulations (Berneche 
and Roux, 2000) have all indicated that alkali ions pre- 
ferentially occupy alternating binding sites within the 
filter, with water molecules interspersed in between. 
This arrangement avoids the strong electrostatic repul- 
sion that would occur when ions occupy neighboring 
sites, an effect that intuitively should be even larger for 
the interaction between the divalent Ba^* ion and an 
alkali lock-in ion. 

To evaluate the proposed scenario within the frame- 
work of molecular dynamics (MD) simulations, we cal- 
culated the PMF of the lock-in K^ ion moving along the 
Z axis of the filter from site Si to the external solution 
when there is a Ba^^ ion occupying site S? (Fig. 2). For 
comparison, we also show the PMF for this ion when site 
S2 is occupied by a K^ ion. The PMF between 2k = 3 A 
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and Zk = 5 A corresponds to the ion occupying 
site Si- This is a highly unstable arrangement when site 
S2 is occupied by Ba^*, with these configurations lying 
>10 kcal/mol higher on the free energy surface. This 
configuration is far less stable in comparison to the con- 
figuration where sites S2 and S4 are occupied by ions, 
in keeping with the higher electrostatic repulsion the 
lock-in ion experiences with a divalent Ba^^ ion in com- 
parison to a ion. This effect diminishes when the lock- 
in ion is in a more external binding site; the binding 
energy of site So is comparable for the Ba^^ and the 
occupied filters, although the lock-in ion tends to 
bind more outwardly in site So when Ba'* is present. 
These PMFs show that is unlikely to bind in site Si 
when site S2 is occupied by Ba^^. A scenario where site Sg 
is the external Ba^* binding site and Si is the binding 
lock-in site is ruled out based on these considerations. 

PMF calculations of barium permeation 

We considered the alternative scenario where binding 
to site So is responsible for blocking the forward translo- 
cation of Ba^^ from site S2 to Si. Previous computational 
studies found that site So is not selective for because 
it is wider than the lower sites, and ions bound in it are 
partially exposed to the external solution (Noskov et al., 
2004; Noskov and Roux, 2006); however, a computa- 
tional study by Kim and Allen (2011) found that site So 
is selective for when the filter is occupied by two Ba"^ 
ions. To comprehensively evaluate if bound in site So 
can block the permeation of Ba^^, we computed the 2D 
PMF of Ba'* translocating from site S2 to Si when there 
is a ion external to it. To evaluate the ion selectivity 
of this lock-in effect, we also computed the PMF when 
the external ion is Na*. In each case, an unrestrained 




Zk(A) 

Figure 2. PMF of a ion along the Z axis ranging from site Si 
(3 A < <6 A), site So (6 A < Zk < 9 A), to the external solution 
(Zk > 10 A) . The green line corresponds to a configuration where 
sites S4 and S2 are occupied by ions, while the orange line cor- 
responds to a configuration where site S4 is occupied by K* and 
site S2 is occupied by Ba^*. 



ion is present in the bottom of the filter, occupying site 
S4 or S3, depending on the configuration of the other ions. 
The 2D PMFs for the translocation of Ba^"^ from S2 to 

51 with an external lock-in ion are presented in Fig. 3. 
In these 2D free energy maps, the x axis corresponds to 
the position of the Ba^*, moving from site S2 (0 A < Zsa < 
2 A) to site Si (3 A < < 5 A), while the y axis corre- 
sponds to the position of the lock-in ion, (left) or Na* 
(right). All ion positions are defined relative to the cen- 
ter of mass of the filter. The PMF for the translocation 
of Ba^^ from S2 to Si shows three distinct minima corre- 
sponding to metastable configurations: Ba^* bound in 

52 with the lock-in ion in So (Fig. 3, bottom left), the 
Ba^^ bound in S2 with the lock-in ion is the external so- 
lution (Fig. 3, top left), and Ba^^ bound in Si with the 
lock-in ion in the external solution (Fig. 3, top right). 

The physical basis for the lock-in effect is immediately 
apparent from these free energy surfaces; the barriers 
are extremely high for direct translocation when the 
lock-in ion is present in So (Zl > 0 A) . This reflects that 
the translocation of Ba^^ to site Si while occupies So 
leads to an arrangement where the ions occupy neigh- 
boring sites and experience a strongly repulsive electro- 
static interaction. The lowest free energy path for Ba^* 
translocation corresponds to the complete exit of the 
lock-in ion into the external solution followed by the 
translocation of Ba^* from S2 to Si. The lowest barrier 
for this transition relative to the minimum is 17 kcal/mol, 
which is consistent with a slow translocation rate of Ba^^ 
in comparison to alkali ions. After crossing this barrier, 
Ba^^ is bound in site Si and the lower unrestrained ion 
has progressed to site S3. The free energy of the state 
where Ba^* is bound in S2 is within 1 kcal/ mol of the 
state where Ba^^ is bound in Si. This indicates that sites 
Si and S2 have similar affinities for Ba'^ in the multi-ion 
binding scenario. However, once the Ba^* is bound in Si, 
external ions can no longer impede its progress toward 
the exit of the pore, as there are no binding sites be- 
yond So. For this reason, within this model, the lock-in 
effect reflects the critical attempt of Ba^^ to translocate 
from S2 to Si. 

The global energy minimum of these surfaces (near 
Zfla = 3.5 A and Zl > 12 A) corresponds to Ba^* occupy- 
ing site Si and the lock-in ion in the external solution. 
The energy minimum of the state where Ba^* occupies 
site S2 and occupies site So is 7 kcal/mol higher in 
energy. This may be an underestimate of the binding 
affinity of site So, as the nonpolarizable force field mod- 
els used in these studies can underestimate the inter- 
action between ions and carbonyls (Yu et al., 2010). 

The PMFs of the Na* and lock-in ions are very simi- 
lar when the ions are outside of site Sq, which is consis- 
tent with a scenario in which the ion is in the external 
solution and therefore effectively uncoupled from the 
filter. The interesting distinctions occur when the lock- 
in ion is in So ( Zl < 9 A) and Ba'* is in S2 ( Zsa < 2 A) . The 
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minimum on PMF occurs when the ion occupies the 
middle of So, near Zl = 7.5 A. The minimum of the PMF 
for Na^ occurs lower, at 2^^ = 7 A, where it can interact 
closely to the four backbone carbonyls. This portion of 
the PMF is '^2 kcal/mol higher in energy for Na^ than 
it is for K*^, which is consistent with site So acting as a 
selective site. 

A recent computational study by Kim and Allen 
(2011) also calculated the PMF of the and Na* ions in 
positions between site So and the external solution while 
the filter was occupied by Ba'*, and proposed that the 
lock-in phenomenon results from binding in site So 
while Ba'* binds in site Sg. Our results are generally con- 
sistent with theirs; however, it should be noted that their 
simulations differs from ours in that two Ba"* ions were 
included simultaneously in the S4 and Sg sites in their 
simulations. Although Ba^^ ion densities were observed 
in both these sites in the structure obtained after grow- 
ing the crystals in 5 mM BaCla (Lockless et al., 2007), 
Piasta et al. (2011) determined that the channel is oc- 
cupied by only one Ba"* ion under the conditions of 
their blockade experiments (from the K*-free block time 
distribution). Based on these observations, the models 
used in our calculations were constructed in such a way 
that there is only one Ba^^ ion in the filter. 

FEP studies of site selectivity 

A comparison of the computed PMFs with the and 
Na^ lock-in ions indicates that site So is moderately selec- 
tive for over Na^ when Ba^^ is bound in site Sg. This 
is somewhat surprising, as previous computational stud- 
ies found that site So is weakly selective for Na^ when 



the lower sites of the filter were occupied by (Noskov 
and Roux, 2006; Egwolf and Roux, 2010). To resolve 
these differences, we quantified the thermodynamic se- 
lectivity of site So when K^, or Ba"*, occupies Sg. This was 
achieved by calculating AAGK-^Na, which is the free en- 
ergy of substituting a bound ion for a Na^ ion relative 
to this substitution in bulk water. 

The free energy of selectivity of a site, AAGK-^Na, is dif- 
ficult to determine accurately using PMFs, as the resolu- 
tion of the free energy surfaces is limited and sampling 
error can cause sizable deviations. Furthermore, this 
method relies on a stable baseline corresponding to the 
ion in the bulk, which requires extending the PMF for 
the ion reaching a large distance away from the pore. 
Computationally, it is more advantageous and straight- 
forward to calculate the free energy of this exchange 
using alchemical FEP. This method computes the free 
energy of converting from one potential energy func- 
tion to another. In this instance, we are computing the 
free energy difference of replacing the ion in site So 
with a Na* ion. We performed two FEP simulations: one 
where site S2 is occupied by a Ba^^ ion and one where it 
is occupied by a ion. These free energies of selectivity 
differ from those reported in the previous studies by 
Noskov and Roux (2006) because we use flat bottom 
restraints and the simulation temperature used here is 
25°C, which is the temperature used in the blockade 
experiments of Piasta et al. (2011), as opposed to Noskov 
and Roux (2006), where the physiological temperature 
of 37°C was used. The ion filter configurations and the 
corresponding free energy of selectivity are illustrated 
in Fig. 4. The details of these calculations are presented 
in the Materials and methods section. 




Figure 3. Contour plots of 2D 
PMF of the permeation of Ba^* 
from site S2 to site Si in the pres- 
ence of either K* (left) or Na"^ 
(right) external to Ba^"^. Neigh- 
boring contour lines differ by 
1 kcal/mol. Low-energy regions 
are depicted in blue and high- 
energy regions are depicted in 
red. The Zg^ axis corresponds to 
the position of the Ba^* ion along 
the Z axis of the system while the 
•Zr/ Zf^a axes correspond to the 
position of the lock-in ion. The 
Z axis is parallel to the pore 
formed by the selectivity filter with 
its zero defined at the center-of- 
mass of the backbone atoms of 
selectivity filter residues 75 — 78. 
Images generated using the crys- 
tal structure (PDB accession no. 
1K4C) are used to indicate the 
position of the sites within the fil- 
ter and their occupancy by Ba^* 
ions (orange), K* ions (green), 
and Na"^ ions (purple) . 
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Consistent with previous studies, site So is essentially 
nonselective when and S4 are occupied by K^, with a 
AAGx-^Na = —0.5 kcal/mol. Remarkably, this value 
changes significantly when a Ba^^ ion occupies site Sg, 
leading site So to become modestly selective for K*, with 
a AAGK-^Na = +1-8 kcal/mol. This is consistent with the 
2D PMFs presented in the "PMF of lock-in ion binding" 
section, which are '^2 kcal/mol higher for Na* in the 
areas corresponding to the lock-in ion occupying site 
So. Therefore, the presence of Ba'^ is responsible for a 
shift of ^2.3 kcal/mol. 

To understand the origin of this difference, we calcu- 
lated the axial distribution function of the lock-in ions 
when bound in site So (Fig. 5). When occupies site 
S2, Na* occupies a broad range of positions within the 
site, but tends to bind at the inner edge of So, where it 
can coordinate directly to the four backbone carbonyls, 
a feature noted previously (Shrivastava et al., 2002; Kim 
and Allen, 201 1 ) , while K* tends to bind in a more outer 
fashion in the site (further away toward the extracellu- 
lar side) , coordinating with both the Tyr78 carbonyl li- 
gands and water molecules from the external solution. 
When Ba^^ occupies Sg, both ions are pushed outward, 
with most probable positions near Z= 7.5 A. But the im- 
pact on the coordination environment of is minimal. 
In contrast, Na* is shifted to a significantly more outer 
position; its coordination environment becomes more 
similar to that of K^, which causes this site to become 
selective when Ba^^ occupies site S2. The trend from the 




K^Na 



-0.5 ± 0.1 kcal moM 




'K^Na 



+1.8 ± 0.1 kcal moM 



Figure 4. Calculated ion selectivity for the S4 K^/Sg K*/So and 
S4 K"^/S2 Ba^*/So ion configurations. AAGK->Na is the free energy 
of replacing the K* ion with a Na* ion, relative to this substitution 
in water. 



free energy simulations is in qualitative agreement with 
the experimental estimates of Piasta et al. (2011). How- 
ever, there is a considerable quantitative discrepancy, as 
the calculated AAGx-^Na of +1 .8 kcal/ mol is considerably 
smaller than the experimental estimate of >7 kcal/mol. 

PMF interpretation of barium blockades 
The 2D PMFs display features that are qualitatively con- 
sistent with the concept of a Ba^^ blockade that is affected 
by the presence of an external lock-in ion, although this 
description does not directly demonstrate how the rate 
of Ba^* exit into the external solution corresponds to 
these PMFs. It is useful to clarify this matter further by 
calculating the rate of permeation of the Ba^^ ion from 
site S2 to Si, which can be described using transition 
state theory (Glasstone et al., 1941; Chandler, 1978) , 

sl/2 -\W*/k^T 



2nm 



-[W(Z')-W(Z„)]/*bT 



(6) 



where Wis the PMF corresponding to the translocation 
of the Ba"^ outwards from the filter along reaction coor- 
dinate Z, AVF* is the height of the PMF relative to the 
reactant well, Z„ is the location of the minimum corre- 
sponding to Ba^^ occupying site S2, Z* is the location of 
the maximum, and k is the transmission coefficient (0 < 
K < 1) . The calculation of these terms is described in the 
Materials and methods section. 

Although Eq. 1 is a familiar expression, the effect of the 
lock-in ion is implicitly included in the activation free en- 
ergy AW*. This is where the difficulty lies because con- 
structing the generalized ID PMF for the permeation of 
Ba^^ is more complicated, as the PMF depends on the 



0.20 




Figure 5. The axial distribution function of Na* and bound in 
site Sfl. The solid lines correspond to the scenario where a ion 
occupies site Sj while the broken lines correspond to Ba^* occupy- 
ing site S2. The mean axial position of the backbone carbonyl oxy- 
gen from residues TyrYS (078) and Gly79 (079) is shown with 
vertical green lines. 
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concentration of the external cation. As the 2D PMFs il- 
Itistrated, the presence of the lock-in ion can consider- 
ably affect the ID PMFs corresponding to Ba^* permeation. 
The probability of an ion being bound to the filter de- 
pends on the external concentration of the lock-in ion, 
so application of transition state theory requires that we 
define a ID PMF of the 2^;, coordinate that varies as a 
function of concentration of the lock-in ion in the exter- 
nal solution. To this end, we express the mean force as a 
linear combination of the mean force when there is no 
external lock-in ion, W„=o(2«a) and the PMF when there is 
an external lock-in ion VKi=i(2Ba). These mean forces are 
weighted by the probability (P) of the lock-in ion being 
present when the Ba^* is present at point Zga- 

" ^^X'^B^. «=0) ^"=0 (^Ba ) + {^){Zs^, n=l) ^ri=\ (^Ba )" 

The probabilities can be expressed in terms of the 
concentration of the lock-in ion in the external solution 
([C]) and the equilibrium constant of binding to the 
lock-in site when Ba^^ is at Zsa (J^q(-^a) ). 

IF) =(F) ^ + (F) 

The derivation of this equation is presented in the Ma- 
terials and methods section. 

In principle, a PMF along the Z axis could be deter- 
mined by integrating Eq. 2, although for our purposes, 
we only need to consider the two extreme cases: (1) the 
[C] = 0 case, when there are no blocking ions in the 
external solution and the only contribution is from 
{F)i^ZY^,n=oy 3^nd (2) when [C] is at a saturating concen- 
trations and (-F)(Zg^,n=i) will dominate. By calculating the 
PMF for these two limiting cases, we can describe the 
lock-in effect. 

We calculated the ID PMFs of Ba~^ permeation sub- 
ject with lock-in ions Na* or external to it from the 
2D PMF of Ba'^ The ID PMF of Ba^* permeation when 
the lock-in ion is in site So, W„=i(ZBa), can be estimated 
by integrating over the regions of the 2D PMF that cor- 
respond to the lock-in ion being present in site So (Zl = 
[5A, 9A]), 

5A 

The PMF when there is no lock-in ion present can be 
estimated by averaging over the regions of the 2D PMF 
where the lock-in ion is too distant to affect permeation 
(Zl= [9 a, 14A]), 

14A 
9A 



We averaged the Na* and PMF to remove spurious 
differences due to sampling error. The ID PMF for this 
unblocked permeation is presented in Fig. 6. When site 
So is not occupied by an ion, permeation of Ba^^ occurs 
with a barrier of roughly 15 kcal/mol. The Eyring reac- 
tion rate calculated using Eq. 6 for Ba^* permeation 
when there is no lock-in ion present is 208 s \ and the 
Grote-Hynes transmission coefficient is k = 0.25, yield- 
ing a rate constant of '^52 s^^, which is in reasonable 
agreement with the experimentally determined rate 
constant of 204 s ^ This is consistent with a slow but 
nonzero rate of permeation of Ba^^ in the absence of a 
lock-in ion. 

The second limiting case, where ion concentration in 
the external solution is sufficiendy high to saturate site 
So, is plotted in Fig. 6 in purple for Na* and green for 
K*. There are minor differences between the Na^ and 
curves; however, the more significant feature is that 
both PMFs are effectively infinite except when Ba'* is in 
its initial state in site S2. This reflects the high repulsion 
between the two ions when they are placed in neighbor- 
ing sites. In this regime, the rate of outward permeation 
will be zero for both types of lock-in ions. This corre- 
sponds to the kinetic model of Piasta et al. (2011), 
which determined that the rate of outward permeation 
of Ba'* is effectively zero when the external concentra- 
tion of was >1 mM. 

Within this framework, the strong lock-in effect ob- 
served in the experiments of Miller and coworkers 
(Neyton and Miller, 1988; Piasta et al., 2011) would 
occur if the equilibrium constant for an ion binding 
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Figure 6. Reduced PMFs for Ba^* ion permeation from site Sj to 
site Si. The black line corresponds to the PMF when the lock-in 
ion is in the external solution, so the lock-in effect is effectively 
zero. The green and purple lines correspond to the PMF when 
either or Na* is bound in site So, respectively. 
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to the lock-in site were large. A selective lock-in 
effect would occur if the equilibrium constant were 
larger for than Na^. Our computations are not quan- 
titatively consistent with the high selectivity and micro- 
molar affinity for that was experimentally observed 
by Piasta et al. (201 1) . The absolute binding affinity of 
K* to site So is is too low according to the calculated 
PMF (Fig. 2) and and the site is only weakly selective 
for over Na* when site Sg is occupied by Ba'^ accord- 
ing to FEP (Fig. 5). 

Possible explanations for the lack of 
quantitative agreement 

Although the present simulations capture some im- 
portant features of a selective lock-in effect of a Ba'^ 
blockade, they are not quantitatively consistent with the 
experiment. Such inconsistencies can be attributed to 
the limitations of various aspects of the computational 
treatment. For instance, as we used a nonpolarizable 
force field, the induced electronic polarization of the 
environment by the ions is neglected in these models. 
This is of particular concern when considering the ef- 
fect of the medium screening the interaction between a 
divalent Ba^^ ion and the lock-in ion. The development 
of a more accurate force field accounting explicitly for 
the effect of induced polarization will, hopefully, allow 
a more satisfactory representation of ion permeation 
(Lamoureux et al., 2003; Lopes et al., 2007; Yu et al., 
2010; Rowley and Roux, 2012). Furthermore, we used 
the high-resolution crystallographic structure of the 
KcsA channel (PDB accession no. 1K4C), which corre- 
sponds to a functional state with a closed intracellular 
gate. Although the processes of interest are taking place 
in the selectivity filter away from the intracellular gate, a 
high-resolution x-ray structure of the channel in the 
open-conductive state might provide a more realistic 
system for the simulations. Lastly, selectivity is governed 
by small free energy differences, and achieving statisti- 
cally converged results is challenging, though feasible 
through enhanced sampling methods. Although we can 
attribute some of the discrepancies between our simu- 
lations to these technical limitations, current explanations 
for the lock-in effect might also need to be examined. 
In particular, the calculated PMFs show that stability of 
Ba^* binding in site Si is comparable to its binding in 
site Sg. 

The possibility that Ba"* might bind to sites in the 
selectivity filter other than sites S<z and S4 will require 
further consideration. Electronic density from the Ba'^ 
ions was observed only for sites S2 and S4 in the crystal- 
lographic structure, whereas no density was observed 
for sites Si and S3 (Lockless et al., 2007). The resolu- 
tion of the Ba'^-occupied x-ray structure is too low to 
unambiguously determine the occupancy of the bind- 
ing sites, so there are several plausible binding con- 
figurations that ought to be considered in analyzing 



the functional data. One interpretation, adopted by 
Piasta et al. (201 1) , is that Ba^^ does not bind favorably 
in sites Si and S3. Assuming that '^10% occupancy is 
a reasonable minimum threshold for detection of ion 
binding, this would imply that the energy of Ba^^ in 
sites Si and S3 is less favorable by at least 1 .4 kcal/ mol 
relative to sites Sg or S4. Thus, a very small energy differ- 
ence could explain the binding pattern observed in the 
x-ray structure. Alternatively, the observed electronic 
density might reflect a structure where two Ba'* ions 
are bound simultaneously in sites Sg and S4. The x-ray 
structure was obtained from crystals that were grown 
from a solution with 5 mM BaCL, which is considerably 
larger than the micromolar Ba^^ concentrations used in 
the single-Ba^^ occupancy blockade experiments (Piasta 
et al., 2011). The implication is that single occupancy 
of site S3 may be energetically accessible, though not 
observed in the x-ray structure due to the high electro- 
static penalty associated with divalent ions occupying 
adjacent sites. Greater clarity about the occupancy of 
the Ba^^ in the binding sites is essential to definitively 
interpreting these experiments. 

Conclusion 

The simulations in this study were undertaken to help 
interpret the Ba^^ block experiments of the KcsA chan- 
nel performed by Piasta et al. (2011). These calcula- 
tions were designed to examine the translocation of the 
Ba^* ion from site S2 to site Si and how this specific pro- 
cess is affected by the binding of extracellular cations to 
the filter. This choice was motivated by two observations. 
First, the Ba^^ ions in the Ba^^-soaked crystal structure of 
KcsA bind either to site S2 or site S4 of the selectivity fil- 
ter Second, the electrophysiological data suggested that 
the dominant situation during the prolonged blocks oc- 
curs when Ba^^ is bound to the outermost of these two 
sites. Under these conditions, the external lock-in site 
must either be site Si or site So, two possibilities that were 
investigated computationally. 

Using umbrella sampling PMF calculations, we deter- 
mined the affinity of K^ to sites So and Si when site S9 is 
occupied by either Ba^^ or K^. The computations indi- 
cate that binding of a cation (K* or Na*) in site Si while 
a Ba^^ is bound in site S2 is energetically prohibitive. 
Although site Si has a weak affinity for K* when S2 is 
occupied by K*, it is highly unfavorable for binding K* 
when S2 is occupied by Ba^^. This suggests that the pro- 
posal of Piasta et al. (2011) that the lock-in effect is due 
to a block of the translocation of Ba^* from site S2 to site Si 
caused by the binding of a K^ ion in the adjacent site Si 
is unlikely. 

If Si is not the lock-in site, the only remaining exter- 
nal site that could serve in this capacity is So. This con- 
clusion is somewhat unexpected and counterintuitive. 
Previous computational studies have indicated that So 
is not selective for K^. In fact, this observation partly 
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motivated Piasta et al. (2011) to rule out the possibility 
that So may be able to act as a lock-in site for the selec- 
tive binding of external cations. However, the PMF 
calculations presented in this paper indicate that So 
becomes considerably more selective when S9 is occu- 
pied by Ba^^ than had been observed in previous studies 
when the filter was occupied by K*. The trend was con- 
firmed by FEP calculations indicating that site So sites 
became significantly more selective for when Ba^^ is 
bound inside the filter, although it is worth noting that 
the computed selectivities are smaller than the values 
determined by the experiments of Piasta et al. (2011). 
We attribute the increase in selectivity to the greater 
electrostatic repulsion between the divalent Ba^^ ion 
and the lock-in ion. This repulsion drives the lock-in ion 
bound to site So in a more outward position, which dis- 
proportionately affects Na^, as it tends to bind favorably 
in a more inward position in the site. This suggests that 
the selectivity derived from Ba^^ blockade experiments 
differs from the selectivity that occurs under physiologi- 
cal conditions. 

To clarify the mechanism of block by the external 
lock-in site, we calculated the rate of Ba^^ permeation 
using Grote-Hynes theory and our calculated ID PMF 
of Ba^^ translocation from S9 to Si. The calculated rate 
constant for this process is ^52 s \ which is in reason- 
able agreement with the experimental rate of 204 s ^ A 
statistical mechanical formulation of the rate was elabo- 
rated to incorporate the concentration dependence of 
the external lock-in ion. When there is a saturating con- 
centration of the lock-in ion in the external solution, 
the PMF for permeation has an extremely high barrier. 
This reflects a scenario where Ba^* permeates from S9 to 
Si when there is an alkali ion in So. The weighting of 
these PMFs depends on the equilibrium constant for 
the binding of the lock-in ion from the external solu- 
tion, where a high binding affinity will increase the 
weighting of the "locked-in" PMF. Nevertheless, although 
the calculated Ba^* permeation PMFs and estimated 
rate constant were consistent with this kind of lock-in 
scenario, the affinity of to site So and the small selec- 
tivity of over Na* in this site are not in quantitative 
agreement with the experiment. 

The lack of quantitative agreement between the cal- 
culations and the experimental data is not surprising 
and confirms that computational models are imperfect 
and need to be improved. Various limitations and de- 
ficiencies, such as sampling error, inaccuracy of non- 
polarizable force fields, and the appropriateness of the 
crystallographic structures used in the simulations and 
how they correspond to the functional data should also 
be considered. Although the present study was limited 
to an examination of the translocation of the Ba^^ ion 
from site Sg to site Si, it is possible that the observed 
lock-in effects actually involve translocation steps with 
Ba^^ bound to sites that are more inward (e.g., S3 or S4) . 



Ultimately, this will require a complete characterization 
of each translocation step involved when a Ba'^ ion en- 
ters the filter from the intracellular side and proceeds all 
the way across the pore to exit on the extracellular side. 
In our view, such an ambitious undertaking would be 
premature at this point given the current limitations of 
the force field and computational model. Gradual prog- 
ress is being made in these areas, which provides the 
prospect that we will eventually be able to achieve quan- 
titative agreement between this experimental electro- 
physiological data and molecular dynamics simulations. 
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